#include "libsec4.h"
#include "stdio.h"
#include "math.h"
#include "stdlib.h"

static double b, g;
void duffing(double x, double y[MMAX], double dy[MMAX])
{
    dy[0] = y[1];
    dy[1] = b * cos(x) - g * y[1] - y[0] * y[0] * y[0];
}

void v2(double x, double y[MMAX], double dy[MMAX])
{
    dy[0] = y[1];
    dy[1] = -4 * PI * PI * y[0];
}

int main()
{
    printf("================求解Dufing模型==============\n");
    {
        g = 0.1;
        b = 0.5;
        int n = 1000;
        int m = 2;
        double h = 0.1;
        double y0[2] = {0, 1};
        double *y = (double *)malloc(n * m * sizeof(double));
        RungeKutta(duffing, y0, y, m, n, h);

#ifdef SAVEFILE
        double *plot;
        int m2 = 2;
        plot = (double *)malloc(n * m2 * sizeof(double));
        for (int i = 0; i < n; i++)
        {
            *(plot + i * m2) = y[i * m];
            *(plot + i * m2 + 1) = y[i * m + 1];
        }
        char filename[256] = "exercise_tmp1.txt";
        printf("求duffing方程的结果保存在文件%s中 \n", filename);
        savef(filename, plot, n, m2, m2);
        free(plot);
#endif
    }

    printf("================求解两种方式求解初值问题==============\n");
    {

        int n = 1000;
        int m = 2;
        double h = 0.005;
        double y0[2] = {1, 0};
        double *y = (double *)malloc(n * m * sizeof(double));
        double *plot;
        int m2 = 3;
        plot = (double *)malloc(n * m2 * sizeof(double));
        RungeKutta(v2, y0, y, m, n, h); //RungeKutta计算
        for (int i = 0; i < n; i++)
        {
            *(plot + i * m2) = i * h;
            *(plot + i * m2 + 1) = y[i * m];
        }

        double s[n];
        double q[n];
        for (int i = 0; i < n; i++)
        {
            s[i] = 0;
            q[i] = 4 * PI * PI;
        }
        y0[1]=1;
        Numerov(q, s, y0, y, n, h); //Numerov计算
        for (int i = 0; i < n; i++)
        {
            *(plot + i * m2 + 2) = y[i];
        }
#ifdef SAVEFILE
        char filename[256] = "exercise_tmp2.txt";
        printf("求duffing方程的结果保存在文件%s中 \n", filename);
        savef(filename, plot, n, m2, m2);
        free(plot);
#endif
    }
}